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The spreading of an infectious disease can trigger human behavior responses to the disease, which in turn plays 
a crucial role on the spreading of epidemic. In this study, to illustrate the impacts of the human behavioral 
responses, a new class of individuals, , is introduced to the classical susceptible-infected-recovered {SIR) 
model. In the model, state represents that susceptible individuals who take self-initiate protective measures 
to lower the probability of being infected, and a susceptible individual may go to state with a response rate 
when contacting an infectious neighbor. Via the percolation method, the theoretical formulas for the epidemic 
threshold as well as the prevalence of epidemic are derived. Our finding indicates that, with the increasing 
of the response rate, the epidemic threshold is enhanced and the prevalence of epidemic is reduced. The 
analytical results are also verified by the numerical simulations. In addition, we demonstrate that, because 
the mean field method neglects the dynamic correlations, a wrong result based on the mean field method 
is obtained—the epidemic threshold is not related to the response rate, i.e., the additional state has no 
impact on the epidemic threshold. 
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Human behavioral responses to the spreading of 
epidemics have been recognized to have great in¬ 
fluence on the epidemic dynamics. Therefore, it 
is very important to incorporate human behaviors 
into epidemiological models, which could improve 
models’ utility in reflecting the reality and study¬ 
ing corresponding controlling measures. How¬ 
ever, analytically-grounded approaches to the 
problem of interacting effect between epidemic 
dynamics and human behavioral responses are 
still lacking so far, and what can help us bet¬ 
ter understand the impacts of human behavioral 
responses. In this work, by introducing a new 
state into classical SIR model to mimic the 
situation that, when susceptible individuals are 
aware of the risk of infection, they may take self- 
protective measures to lower the probability of 
being infected. We then derive theoretical formu¬ 
las based on the percolation method for the epi¬ 
demic threshold and the prevalence of epidemic. 
Our results indicate that the introduction of the 
new S^ state can significantly enhance the epi¬ 
demic threshold and reduce the prevalence of epi- 
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demic. It is worth mentioning that a wrong re¬ 
sult may be obtained when using the traditional 
mean field method the epidemic threshold is not 
altered by the additional S^ state. The result 
highlights that if the effects of human behav¬ 
ioral responses are ignored in mathematical mod¬ 
eling, the obtained results cannot really reflect 
the spreading mechanism of epidemics among hu¬ 
man population. 


I. INTRODUCTION 

Since network models can well describe the spread¬ 
ing of infectious disease among populations, many efforts 
have been devoted to studying this fieldii^. At first, re¬ 
searchers mainly paid attention to analyze the impact 
of the network structure on epidemic spreading and the 
control strategies, for example, how the network topol¬ 
ogy affects the epidemic threshold and the prevalence of 
epidemio^^— , or how to design an effective immunization 
strategy to control the outbreaks of epidemics^ii^. In re¬ 
ality, outbreaks of epidemics can trigger spontaneous be¬ 
havioral responses of individuals to take preventive mea¬ 
sures, which in turn alters the epidemic dynamics and 
affects the disease transmission processii^— . Thus, re¬ 
cently, some studies have made attempts to evaluate the 
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impact of the human behaviors on epidemic dynamics. 
For instance, Funk et studied the impacts of aware¬ 
ness spread on both epidemic threshold and prevalence, 
and they found that, in a well-mixed population, spread 
of awareness can reduce the prevalence of epidemic but 
does not tend to affect the epidemic threshold, yet the 
epidemic threshold is altered when considering on social 
networks; Sahneh et al. considered a Susceptible-Alter- 
Infected-Susceptible {SAIS) modeli^, and they found 
that the way of behavioral response can enhance the 
epidemic threshold; Melon! et al. constructed a meta¬ 
population model incorporating several scenarios of self- 
initiated behavioral changes into the mobility patterns of 
individuals, and they found that such behavioral changes 
do not alter the epidemic threshold, but may increase 
rather than decrease the prevalence of epidemio^S. Mean¬ 
while, in Ref.—, by designing the transmission rate of epi¬ 
demic as a function of the local infected density or the 
global infected density, Wu et al. investigated the effect 
of such a behavioral response on the epidemic threshold. 

One fact is that, the infectious neighbors can infect a 
susceptible individual, they can also trigger the aware¬ 
ness of the susceptible individua l^^’^^ . In view of this, in 
Ref.—, Perra et al. introduced a new class of individuals, 
, that represents susceptible people who self-initiate 
behavioral changes that lead to a reduction in the trans- 
missibility of the infectious disease, into the classical SIR 
model, and they found that such a model (SS^IR) can 
induce a rich phase space with multiple epidemic peaks 
and tipping points. However, the network structure was 
not incorporated into these models. As we know, when 
the model is considered within the network based frame¬ 
work, new theoretical tools should be used and new phe¬ 
nomena may be observed. In view of this, we incorpo¬ 
rate the network structure into the SS^IR model^^ to 
investigate its spreading dynamics. In the model, when 
contacting an infectious neighbor, susceptible individuals 
may be infected (from S state / state) with a transmis¬ 
sion rate or a behavioral response may be triggered (from 
S state to state) with a response rate. We provide a 
theoretical formula for the epidemic threshold as well as 
the prevalence of epidemic via the percolation method^, 
our results show that the introduction of class can en¬ 
hance the epidemic threshold and reduce the prevalence 
of epidemic. We also demonstrate that a wrong result 
can be obtained—the introduction of class cannot al¬ 

ter the epidemic threshold when using mean field method 
to such a model. The reasons are presented in Sec. El 


II. DESCRIPTIONS OF THE MODEL 

For the classical SIR model on complex network, 
where each node on network can be in one of three states: 
Susceptible {S), Infected (/) or Recovered (R). The 
transmission rate along each SI link is (3, and an infected 
node can enter R state with a recovery rate /i. To reflect 
the fact that, upon observation of infection, susceptible 


individuals may adopt protective measures to lower their 
infection risk, a new class, denoted by S^, is introduced 
into the original SIR model, we use SS^IR model to 
denote the modified SIR model in this study. In the 
model, when an S node contacts an / neighbor, besides 
the probability of being infected, the S node can go to 
S^ state with a response rate I3p > 0. The transmission 
rate for the S^ nodes is lower, thus we assume the trans¬ 
mission rate for S^ nodes is 7 / 3 , where 0 < 7 < 1 is a 
discount factor. 

The SS^IR model is described by the four following 
reactions and the associated rates: 

5 + /A 2/, (1) 


S + I S^ + 1, (2) 


S^ + 1^ 2/, (3) 


/ A R. (4) 

Note that the SS^IR model returns to the SIR model 
once Pf = 0 , and the S^ state corresponds to fully vac¬ 
cinated state when 7 = 0 . 

III. THEORETICAL ANALYSIS 

In our model, during a sufficiently small time interval 
At, the transitionrates of an SI edge becoming an //, 
S^I and SR edge are /3, /3f and /r, respectively. As a 
result, the probabilities of an SI edge becoming an II and 
S^ I edge are given as A - p+Pp+f and r 2 - 
respectively. Similarly, since the transition rate of an 
S^I edge becoming an II and S^R during a sufficiently 
small time interval At are 7/3 and /i, the probability of 
an S^I edge becoming an II edge is 

To analyze our proposed model, we first de¬ 
fine “externally infected neighbor” (EIN) for any node. 
For node i, if a neighbor j is an EIN means that j is 
infected by its neighbors other than i ( i.e., j is infected 
even though node i is deleted from the networks, which 
is the basic assumption of the cavity theory in statisti¬ 
cal physics. Note that this method is suitable for the 
networks with negligible number of loops as the network 
size is sufficiently large^^). The probability of neighbor 
j being an EIN of i is defined as 9, then, for the a node 
with degree fc, the probability of having m EINs is given 
as: 

p{m\k) = Qe^{i-ef-'^ 


( 5 ) 
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Let p{R\m) be the probability of i being infected when 
it has number m of EIN. To calculate such a probability, 
we need to recognize that, in our model, an S node can 
become an I node through two ways: (a) the S node 
is directly infected; or (b) the S node first goes to 
state and then is infected. To facilitate the analysis, we 
approximately assume that the impacts of Ts infected 
neighbors on node i happen in a non-overlapping order, 
i.e., they play roles on node i one by one. 

For case (a), the probability of node i being infected 
by the sth infected neighbor is given as: 


with degree k being infected is 

p{R\k) = p{R\m)p{m\k) = 1 — (1 — 9Ti — 6 T 2 ) 

m 

(1 - T^)T 2 


= 1 - 


Ti + T 2 - Ta 
Ti - T 3 + r2T3 


T1+T2- T3 

(1 - n)T 2 

'T 1 + T 2 -T 3 


[(1 - eTs)'^ - (1 - 6T3 - 0x2)^] 
(1 - 0Ti - 0T2f 

(1 - er^f. (12) 


= (i-ri-T2)*-iTi, 


Then the EIN probability 9 is the solution to the self- 
consistent condition 


Eq. ([ 6 ]) indicates that previous s — 1 infected neighbors 
have not changed the state of node i (not become I or 
state) before they become R state. Therefore, the 
probability of i being infected is: 

m rp 

P3iR\m) = ^ - (1 - Ti - r2)-]. (7) 

For case (b), node i should first become state, and 
the probability that the susceptible node i is altered by 
the Ith infected neighbors and becomes can be written 
as: 


9 = J 2 Q{k)p{R\k) 


= 1 - 


Ti - T 3 + T 2 T 3 


T 1 +T 2 - T 3 
(1 - T3)T2 
'T 1 +T 2 -T 3 


Gi{l - 9 Ti - 9T2) 
Gi{l- 9 T 3 ) = f{ 9 ). (13) 


In Eq. (IT^ . Q{k) = jg excess degree dis¬ 

tribution, where P{k) is the degree distribution and (fc) 
is the average degree. The generating function for Q{k) 
is given as 


7 I 2 = (1 - Ti - T2)'-1T2, (8) 


Gi(a:) = ^Q(fc)x^ (14) 


which also indicates that previous I — 1 infected neighbors 
have not changed the state of node i before they become 
R state. For the remainder m — I + 1 infected neighbors 
(including the infected neighbor who just made i go to 
state), they can infect i with probability: 

A 3 = i-(i-r3)™-'+i. (9) 

As a result, the probability of node i first becoming 
state and then going to I state is: 

m rp 

P 2 {R\m)= Y^A2*A3 = [1 - (1 - Ti - T 2 )n 

i=i 1+2 

- tI+t2^-t3 ^^^ ~ - (1 - - T2r]-m 


The probability p{R\m) is 
?|r? 

[(1 - T 3 T - (1 - Ti - T2r] 


p{R\m) = pi{R\m) -f p 2 {R\m) = l-{l-Ti-T 2 ) 
(1 - T3)T2 


= 1 - 


T 1 +T 2 - T 3 
T 1 -T 3 + T 2 T 3 


T1+T2- T3 
(1 - T3)T2 
'T3+T2-T3 


(1 - Ti - T 2 )'" 

(1-T3)™. (11) 


Combing Eqs. m and (EH), the probability of a node 


There is a trivial solution 0 = 0 in self-consistent equa¬ 
tion m- To have a non-trivial solution, the following 
condition must be met: 


^^|,=o = (Ti+r2T3)Gl(l)>l, (15) 

which implies the epidemic can outbreak when 

(3 /r + 7/3 + 7/3 f ^ (fc) 
p + jP p + /3 + Pf ~ (k^) - (k)' 


For the prevalence of epidemic (defined as i?(oo)), 
we can numerically solve 9 from self-consistent equa¬ 
tion EH. then the formula of R{oo) is 


i?(oo) = ^P(fc)p(i?|fc) 

Ti - T 3 + T 2 T 3 


k 

= 1 - 


T 1 +T 2 - T 3 
(1 - T3)T2 
Ti -tT 2 -T 3 


Go(l - 9Ti - 0 X 2 ) 
Go(1-0T3). (17) 


In Eq. (fT71) . Go(.) is the generating function of degree 
distribution P{k), which is described as: 


Go(a:) = ^P(fc)x^ 

A;=0 


(18) 
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IV. SIMULATION RESULTS 

In this section, we perform an extensive set of Monte 
Carlo simulations to validate the theoretical predictions 
in Section Um Here we carry out simulations on an Erdos- 
Renyi network (labeled ER network)^ with network size 
N = 10000 and average degree (k) = 10, and a config¬ 
uration network generated by an uncorrelated configu¬ 
ration model (UCM)21. The configuration network also 
has N = 10000 nodes and the degree distribution meets 
P{k) ^ whose minimal and maximal degrees are 

krmn = 3 and kmax = VN , respectively. 

A. Results on ER network 



A P 


Differing from the SIS (Susceptible-Infected- 
Susceptible) model, it is not an easy thing to determine 
the epidemic threshold for the SIR model owing to the 
non-zero value of R. In doing so, in Refi^, Shu et al. 
suggested that the variability measure 


Vip^) - {pV 

ip) 


(19) 


can well predict the epidemic threshold for the SIR 
model, where p denotes the prevalence of epidemic in 
one simulation realization^Si^. A can be explained as 
the standard deviation of the epidemic prevalence, and 
is a standard measure to determine critical point in equi¬ 
librium phase on magnetic systent^. In our simulations, 
we have taken at least 1000 independent realizations to 
predict the epidemic threshold. For convenience, in this 
study, we set recovery rate p = 1 . 0 . 

In Fig. [H for different response rate /3 f, the value of 
i?(oo) (upper panels) and the measure A (lower panels) 
as the functions of the transmission rate /3 are investi¬ 
gated. As shown in Fig. [U one can observe that the 
variability measure can well predict the epidemic thresh¬ 
old for our SS^IR model. As a result, in the following 
figures, we use this method to determine the epidemic 
threshold, i.e., the point where the value of A is the 
maximal. Fig. [T]also describes that, no matter 7 = 0.1 
[see Figs. [TKa)-(b)] or 7 = 0.3 [see Figs. [Il[c)-(d)], on 
the one hand, the epidemic threshold is enhanced as the 
response rate /3 f is increased. On the other hand, for 
the a fixed value of /I, Figs. [T][a)and (c) suggest that the 
prevalence of epidemic is remarkably reduced when /3 f 
is increased. The results suggest that, by introducing 
an additional protective state, S^, to the classical SIR 
model, the conclusions are quite different from the pre¬ 
vious results which have not incorporated the impacts of 
human behavioral responses. The result again empha¬ 
sizes the fact that the spontaneous behavioral responses 
of individuals to the emergent diseases have vital impacts 
on the epidemic dynamics. If the behavioral responses 
are ignored in mathematical modelling, the obtained re¬ 
sults cannot really reflect the spreading mechanism of 
epidemics among human population. 


FIG. 1. (Color online) On ER networks, the prevalence of 
epidemic R{oo) (upper panels) and the variability measure 
A (lower panels) as functions of the transmission rate j3 for 
different values oi Pf and 7 . Note that, here fip is a, rate 
rather probability, as a result, whose value may be larger than 
one. (a)-(b) 7 = 0.1; (c)-(d) 7 = 0.3. The pink lines are given 
to demonstrate that the peak value of A corresponds to the 
epidemic threshold Pc- 


We then compare the theoretical results with the 
Monte Carlo simulations on ER network in Fig. [2] and 
Fig. [3l Since the degree distribution of an Erdds-Renyi 
network is P{k) = the generating functions 

meet: 


Go(a:) = C?i(a;) = (20) 


According to Ineq. dTSl) . the epidemic threshold Pc for 
ER network is determined by the following equation 

Pc P + iPc + IPf _ 1 . , 

p + 7/^0 p Pc P Pf {k) 

Moreover, substituting Eq. (1^ into Eq. (IT^ and 
Eq. (fT71) . the prevalence of epidemic R{oo) can be easily 
solved. 

The comparison of R{oo) between the simulations and 
the theoretical result is plotted in Fig.[2J which indicates 
that the numerical simulation and the theoretical result 
are in good agreement. Meanwhile, the epidemic thresh¬ 
old for Pc obtained from Eq. (ED) and from numerical 
method (i.e., the point where A is maximal) is compared 
in Fig. [31 which also indicates that the epidemic thresh¬ 
old predicated by our method is remarkable agreement 
with numerical simulations. The result in Fig. [3] also 
suggests that the epidemic threshold Pc is increased as 
the value of 7 is decreased. Importantly, the reduction is 
more efficient when the response rate Pf is larger. 
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FIG. 2. (Color online) Comparison between the Monte Carlo 
based simulations and the theoretical predictions for R{oo) on 
ER networks. The simulation results are denoted by symbols 
and the theoretical predictions are denoted by the correspond¬ 
ing lines. The theoretical results are obtained by substituting 
Eq. (Uni) into Eqs. (HI and dUl), and the value of 6 is numer¬ 
ically solved from Eq. (na). then R{oo) is got from Eq. (ITTll . 
(a) 7 = 0.1; (b) 7 = 0.3. 


FIG. 4. (Color online) Comparison between the Monte Carlo 
based simulations and the theoretical predictions for R{oo) 
on UCM networks. The simulation results are denoted by 
symbols and the theoretical predictions are denoted by the 
corresponding lines. The theoretical results are obtained by 
substituting a fixed degree distribution P{k) into Eqs. (1181) 
and and then R{oo) can be solved from Eqs. (1181) and 
(flSl) as described in Fig. (a) 7 = 0 . 1 ; (b) 7 = 0.3. 



FIG. 3. (Color online) Comparison between the Monte Carlo 
based simulations and the theoretical predictions for the epi¬ 
demic threshold /3c on ER networks. The theoretical predic¬ 
tions denoted by lines are obtained from Eq. m, and the 
simulation results denoted by symbols are the points where 
the values A are maximal, (a) /3c as a function of the response 
rate /3f for different values of 7 ; (b) /3c as a function of the 
discount factor 7 for different values of /3f. 


B. Results on UCM network 

Real complex networked systems often possess certain 
degree of skewness in their degree distributions, typically 
represented by some scale-free topology. We thus check 
our model on UCM network with degree distribution 
P{k) ^ k~^'^ to illustrate that our theory can general¬ 
ize to the networks with heterogenous degree distribution 
and in the absence of degree-to-degree correlation. 

As shown in Fig.Uand Fig.jSJ one can see that the an¬ 
alytical results are in good agreement with the numerics. 



FIG. 5. (Color online) Comparison between the Monte Carlo 
based simulations and the theoretical predictions for the epi¬ 
demic threshold /3c on UCM networks. The theoretical pre¬ 
dictions denoted by lines are obtained from Ineq. dSJ, and the 
simulation results denoted by symbols are the points where 
the values A are maximal, (a) /3c as a function of the response 
rate /3 f for different values of 7 ; (b) /3c as a function of the 
discount factor 7 for different values of /3 f. 


They also indicate that, since increasing Pf can induce 
more susceptible individuals go to state and reduc¬ 
ing 7 can lower the risk of nodes, as a result, both 
of them can lower the prevalence of epidemic R{oo) and 
increase the epidemic threshold Pc- 


V. MEAN FIELD METHOD FOR THE MODEL 

One possible way to describe our proposed SS^IR 
model is the mean field method, which can be written 
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FIG. 6 . (Color online) Based on the differential Eqs. (1221) - 
the prevalence of epidemic iZ(oo) as a function of /3 is 
presented, (a) 7 = 0.1; (b) 7 = 0.3. 


as: 

-f3kSkO -/SphSkO, ( 22 ) 

^FkSkO - l3-fkS^ 0, (23) 

kQP[Sk+iS^)-h, (24) 

4 , (25) 

where the factor 0(t) = P{k'\k)Ik'{t) represents the 

probability that any given link points to an infected 
node. In absence of any degree correlations, 0(t) = 

Based on the above differential equations, we can ob¬ 
tain that the epidemic threshold jSc = (detailed 
derivation is given in section fVIip . which means that the 
epidemic threshold for our model is not related to the 
response rate Pf or the discount factor 7 , and which 
is the same to the epidemic threshold of classical SIR 
model. The simulation results based on the Eqs. (I22II25I1 
in Fig. ini also indicate that, based on mean field method, 
the epidemic threshold is not altered by different values 
of Pf or 7 . However, our previous simulation results 
based on Monte Carlo method have suggested that the 
epidemic threshold j3c is increased when /3 _f is increased 
or 7 is reduced. That is to say, the conclusion obtained 
by mean-field method is wrong. 

Now let us explain why the mean field method gives a 
wrong result. It is known that, near the epidemic thresh¬ 
old, the fraction of infected nodes (label pi) is very small. 
When using the mean field method, the dynamic corre¬ 
lation is neglected, the probability of S node becoming 
is proportional to 0 {pi) since the average fraction 
of infected nodes among neighborhood equals to 0 {pi). 
Similarly, the probability of node becoming I node is 
also proportional to 0{pi). As a result, the probability 
of ^ / is proportional to 0 {p'j), which leads to 

the effect of the is ignored and the epidemic thresh¬ 
old obtained by the mean field method is not related to 
the value of ,5^ or 7 . In fact, when an S node becomes 


dSk{t) 

dt 

dS[{t) 

dt 

dlkjt) 

dt 

dRk{t) 

dt 


an node there must exist at least one infected nodes 
among the neighborhood of the S node. More impor¬ 
tantly, these infected neighbors may exist for a certain 
time interval, so the probability of S' —?► —)► / is not 

proportional to 0{pj). However, as deduced in Eq. ©, 
the dynamics correlation near the epidemic threshold is 
considered in our above analysis, which can accurately 
predict the epidemic threshold. 


VI. CONCLUSIONS 

To summarize, we have proposed an SS^IR epidemi¬ 
ological model in complex networks, in which the prob¬ 
ability of susceptible individuals becoming state is 
proportional to the number of infected neighbors, to re¬ 
flect the fact that individuals are more likely to take 
protective measures when they find their neighbors are 
infected. By using theoretical analysis as well as nu¬ 
merical simulations, we found that the prevalence of epi¬ 
demic is effectively reduced and the epidemic threshold 
is remarkably increased when the response rate Pf is in¬ 
creased or the discount factor 7 is reduced. Moreover, 
we have demonstrated that the mean field based analysis 
provides a wrong result: the epidemic threshold is not 
related to the response rate Pf or discount factor 7 . The 
reason is that, near the epidemic threshold, the probabil¬ 
ity of S' —>■ > J is a second order infinitesimal since 

the mean field method ignores the dynamic correlation, 
which makes the effect of state to be ignored. 

With the development of technology, information in¬ 
duced awareness or self-protective behaviors can not only 
diffuse through the contact networks where the diseases 
spread but also can fast diffuse through many different 
channels, such as, the word of mouth, news media, on¬ 
line social networks, and so on. In view of this, recent 
well-studied multiplex network theory may be an ideal 
framework to mimic the interplay of information or re¬ 
lated awareness and the epidemic dynamics^^^— . Thus, 
how to generalize our model to multiplex networks and 
provide theoretical analysis to the model is a challenge 
in our further work. 


VII. APPENDIX 

When assuming Sk{0) « 1, then from Eq. (l22l) . we 
have^ 

Skit) = e-M/3+/3^)<^(t)^ (26) 

where pit) = f* 0 (t')dt' = j^J^k kPikpRkit). 
Substituting Eq. (l26|) into Eq. (|23)) . one has: 

pqF 

= PFkQit)e-^^^+^^'>‘^^*^ - p-fkS[ 0(t). (27) 
By using the variation of constants method, S[it) is 
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solved as: 

s^it) = 


Then, 


Pf 


P + Pf- Pi 


g-/37fc0(t)(i _g 


-k{P+PF-P~t)<l>{t)\ 


(28) 


Q{t) = = 1 - 

OjL \ fC) 




-P^k<l>{t) 


(1 - e 


-fc(/3+/3F-/37)0(t)') 


(29) 


By letting ^co = limt^oo ^(t), and with the fact that 


J 


= 1 - -^EfcP(fc)e“'=('^+^'")'^'’ 

(fc) 


-P~tk(t> 


The value poo = 0 is always a solution. In order to have 
a non-zero solution, the condition 


dfipoo) I 

—77-l0=o=O 

^0oo 


(/3 + /3;^)(fc2) 


{k) 


P{k^) 

(fc) 


> 1 


PFiP) 

(k) 


(31) 


should be satisfied, which means that the epidemic 
threshold Pc = 
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